---
title: "replication"
author: "Iseki, Tatsuya and Shinomoto, Sou" 
output: html_notebook
---

##Load and attach packages used for replication
If you have not installed packages bellow, please install before.

```{r}
library(readr)
library(dplyr)
library(tidyr)
library(plyr)
library(estimatr)
library(margins)
library(ggplot2)
library(coefplot)
```

##Preprocessing
Load the data for replication.

```{r}

d <- read_csv("data.csv")
View(d)
```

Exclude satisficer and respondents who fail to complete questionnaire.
```{r}
d <-  d %>%   filter(
               !is.na(d$Q1.1) ||
               !is.na(d$Q2.3) ||
               !is.na(d$Q3.1) || 
               !is.na(d$Q16.1) || 
               !is.na(d$Q11.3) ||                                       
               !is.na(d$Q12.3) ||
               !is.na(d$Q13.3) ||
               !is.na(d$Q14.3) ||
               !is.na(d$Q15.3) ||
               !is.na(d$Q18.1) ||
               d$Q1.1 == "同意する" || 　　　　　　　　　　　　　　　　　
               d$Q2.3 == "とても関心がある,ある程度関心がある" ||        
               d$Q3.1 != "18歳未満" ||                                   
               d$Q3.1 != "70歳以上" ||
               d$Q16.1 == "イラン" ||                                    
               d$Q18.1 == "確認した")                                                                 

d <- filter(d, !is.na(Q18.1)) ##ダメ押し

View(d)
```


Prepare the explanatory variables ("veteran" for retired general officer)
```{r}
d$control    <- ifelse(d$Q11.3 != "NA", 1, 0)
d$control    <- replace_na(d$control, 0)     

d$anonymous  <- ifelse(d$Q12.3 != "NA", 1, 0)
d$anonymous  <- replace_na(d$anonymous, 0)

d$opposition <- ifelse(d$Q13.3 != "NA", 1, 0)
d$opposition <- replace_na(d$opposition, 0)

d$scholars   <- ifelse(d$Q14.3 != "NA", 1, 0)
d$scholars  <- replace_na(d$scholars, 0)

d$veteran    <- ifelse(d$Q15.3 != "NA", 1, 0)
d$veteran    <- replace_na(d$veteran, 0)


```


Ideology for hypothesis 2
```{r}
##イデオロギー
d$ideology <- d$Q6.17_1

##イデオロギー三点尺度
d$ideology_left <- ifelse(d$ideology <= 3, 1, 0)
d$ideology_centre <- ifelse(d$ideology <= 3, 0,
                            ifelse(d$ideology < 8, 1, 0))
d$ideology_right <- ifelse(d$ideology >= 8, 1, 0)

d$ideology3 <- ifelse(d$ideology_left == 1, 1,
                      ifelse(d$ideology_centre == 1, 2, 3))
```


Confidence in institutions for hypothesis 3
```{r}
##Confidence in Constitutional Scholars and Universities

d$trust_scholars_pretreatment <- d$Q6.3_7 #constitutional scholars
d$trust_university_pretreatment <- d$Q6.3_6 #Universities


##Confidence in the Diet
d$trust_parliament_pretreatment <- d$Q6.3_2
d <- filter(d, !is.na(trust_parliament_pretreatment))

##Confidence in the SDF
d$trust_SDF_pretreatment <- d$Q6.3_5
d<- filter(d, !is.na(trust_SDF_pretreatment))

```



outcome variable:Support for dispatch
```{r}

d$dispatch_control <-revalue(d$Q11.3, 
                             c("支持する" = 4,
                               "どちらかといえば支持する" = 3,
                               "どちらかといえば支持しない" = 2,
                               "支持しない" = 1))
d$dispatch_control <- replace_na(as.numeric(d$dispatch_control), 0)

d$dispatch_anonymous <- revalue(d$Q12.3, 
                             c("支持する" = 4,
                               "どちらかといえば支持する" = 3,
                               "どちらかといえば支持しない" = 2,
                               "支持しない" = 1))
d$dispatch_anonymous <- replace_na(as.numeric(d$dispatch_anonymous), 0)


d$dispatch_opposition <- revalue(d$Q13.3, 
                             c("支持する" = 4,
                               "どちらかといえば支持する" = 3,
                               "どちらかといえば支持しない" = 2,
                               "支持しない" = 1))
d$dispatch_opposition <- replace_na(as.numeric(d$dispatch_opposition), 0)

d$dispatch_scholars   <- revalue(d$Q14.3, 
                             c("支持する" = 4,
                               "どちらかといえば支持する" = 3,
                               "どちらかといえば支持しない" = 2,
                               "支持しない" = 1))
d$dispatch_scholars <- replace_na(as.numeric(d$dispatch_scholars), 0)

d$dispatch_veteran   <- revalue(d$Q15.3, 
                             c("支持する" = 4,
                               "どちらかといえば支持する" = 3,
                               "どちらかといえば支持しない" = 2,
                               "支持しない" = 1))
d$dispatch_veteran <- replace_na(as.numeric(d$dispatch_veteran), 0)


d$outcome_dispatch <- d$dispatch_control + d$dispatch_anonymous + d$dispatch_opposition + d$dispatch_scholars + d$dispatch_veteran

```


Control variables：SES, demographic traits, etc
```{r}
#Age
d$Dem_Age <- as.numeric(d$Q3.1 )

#Gender（Base category: male）
d$Dem_Gender_Female <- ifelse(d$Q5.1 == "女性", 1, 0)
d$Dem_Gender_Others <- ifelse(d$Q5.1 == "その他", 1, 0)

#Resident prefecture
d$Dem_Prefecture <- as.factor(d$Q5.3)

#Education background
d$SES_Education <- d$Q5.7

d$SES_Education_HighSchool <- ifelse(d$SES_Education == "高校", 1, 0)
d$SES_Education_VocationalHighSchool <- ifelse(d$SES_Education == "短大／高専／専門学校", 1, 0)
d$SES_Education_University <- ifelse(d$SES_Education == "大学", 1, 0)
d$SES_Education_GraduateSchool <- ifelse(d$SES_Education == "大学院", 1, 0)

#Household income
d <- filter(d, d$Q5.9 != "答えたくない")
d$SES_Income <- mapvalues(d$Q5.9, c("300万円未満",
                                  "300~350万円未満",
                                  "350~400万円未満",
                                  "400~450万円未満",
                                  "450~500万円未満",
                                  "500~550万円未満",
                                  "550~600万円未満",
                                  "600~650万円未満",
                                  "650~700万円未満",
                                  "700~750万円未満",
                                  "750~800万円未満",
                                  "800~850万円未満",
                                  "850~900万円未満",
                                  "900~950万円未満",
                                  "950~1000万円未満",
                                  "1000万円以上"
                                  ),
                                  c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16))
d$SES_Income <- as.numeric(d$SES_Income)

d$SES_Income_300under <- ifelse(d$SES_Income == 1, 1, 0)
d$SES_Income_300_350 <- ifelse(d$SES_Income == 2, 1, 0)
d$SES_Income_350_400 <- ifelse(d$SES_Income == 3, 1, 0)
d$SES_Income_400_450 <- ifelse(d$SES_Income == 4, 1, 0)
d$SES_Income_450_500 <- ifelse(d$SES_Income == 5, 1, 0)
d$SES_Income_500_550 <- ifelse(d$SES_Income == 6, 1, 0)
d$SES_Income_550_600 <- ifelse(d$SES_Income == 7, 1, 0)
d$SES_Income_600_650 <- ifelse(d$SES_Income == 8, 1, 0)
d$SES_Income_650_700 <- ifelse(d$SES_Income == 9, 1, 0)
d$SES_Income_700_750 <- ifelse(d$SES_Income == 10, 1, 0)
d$SES_Income_750_800 <- ifelse(d$SES_Income == 11, 1, 0)
d$SES_Income_800_850 <- ifelse(d$SES_Income == 12, 1, 0)
d$SES_Income_850_900 <- ifelse(d$SES_Income == 13, 1, 0)
d$SES_Income_900_950 <- ifelse(d$SES_Income == 14, 1, 0)
d$SES_Income_950_1000 <- ifelse(d$SES_Income == 15, 1, 0)
d$SES_Income_1000over <- ifelse(d$SES_Income == 16, 1, 0)


#Interest in Politics
d <- filter(d, !is.na(Q6.1))
d$SES_Interest <- mapvalues(d$Q6.1, c("まったく関心がない", "あまり関心がない", "ある程度関心がある", "とても関心がある"), c(1,2,3,4))
d$SES_Interest <- as.numeric(d$SES_Interest)



```




##Analyses

#Simple treatment effect: hypothesis 1
```{r}


model1 <- lm_robust(outcome_dispatch ~ anonymous + scholars + opposition + veteran +
                      Dem_Age + Dem_Gender_Others + Dem_Gender_Female +
                      as.factor(SES_Income) + SES_Education + SES_Interest,
                    se_type = "stata", data = d)

summary(model1)

#m <- margins(model1, 
#             variables = c("anonymous", "scholars", "opposition", "veteran"),
#             )

#plot(m)

coefplot(model1, 
         lwdOuter = 1, 
         coefficients = c("veteran",  "opposition","scholars", "anonymous"),
         ylab = "",
         xlab = "Average Treatment Effect",
         newNames = c(annonymous = "Anonymous Criticism", scholars = "Constitutional Scholars", opposition = "Opposition Party", veteran = "Retired General Officers")
         ) +
  theme(text = element_text(size = 16))

```



#Conditional Effect (ideology): hypothesis 2
```{r}

model2.1 <- lm_robust(outcome_dispatch ~ anonymous*ideology + scholars*ideology + opposition*ideology + veteran*ideology,
                      se_type = "stata", data = d) 

model2.2 <- lm_robust(outcome_dispatch ~ anonymous*ideology + scholars*ideology + opposition*ideology + veteran*ideology +
                        Dem_Age + Dem_Gender_Others + Dem_Gender_Female +
                        as.factor(SES_Income) + SES_Education + SES_Interest,
                      se_type = "stata", data = d)

#cplot(model2.1, x = "ideology", dx = "scholars", what = "effect")

cplot(model2.2, x = "ideology", dx = "scholars", what = "effect", rug = FALSE, ylim = c(-0.65, 0.55), 
      main = "Criticism from Constitutional Scholars", xlab = "Ideology", ylab = "Marginal Effect",cex.main =1.5, cex.lab =1.5, cex.lab = 1.5, draw = TRUE) 
abline(h=0, lty=2)

cplot(model2.2, x = "ideology", dx = "anonymous", what = "effect", rug = FALSE, ylim = c(-0.65, 0.55),
      main = "Anonymous Criticism", xlab = "Ideology", ylab = "Marginal Effect", cex.main =1.5, cex.lab =1.5, cex.lab = 1.5, draw = TRUE)
abline(h=0, lty=2)

cplot(model2.2, x = "ideology", dx = "opposition", what = "effect", rug = FALSE, ylim = c(-0.65, 0.55),
      main = "Criticism from Opposition Party",xlab = "Ideology", ylab = "Marginal Effect",cex.main =1.5, cex.lab =1.5, cex.lab = 1.5, draw = TRUE)
abline(h=0, lty=2)

cplot(model2.2, x = "ideology", dx = "veteran", what = "effect", rug = FALSE, ylim = c(-0.65, 0.55),
      main = "Criticism from Retired General Officers",xlab = "Ideology", ylab = "Marginal Effect",cex.main =1.5, cex.lab =1.5, cex.lab = 1.5, draw = TRUE)
abline(h=0, lty=2)
```


#Conditional effect (confidence): hypothesis 3
```{r}


#model3.1 <- lm_robust(outcome_dispatch ~ anonymous*Technocracy_attitude + scholars*Technocracy_attitude + 
#                        opposition*Technocracy_attitude + veteran*Technocracy_attitude,
#                      se_type = "stata", data = d) 

model3.2 <- lm_robust(outcome_dispatch ~ anonymous + 
                        scholars*trust_scholars_pretreatment + 
                        opposition + 
                        veteran +
                        Dem_Age + Dem_Gender_Others + Dem_Gender_Female +
                        SES_Income + SES_Education + SES_Interest,
                      se_type = "stata", data = d)

model3.3 <- lm_robust(outcome_dispatch ~ anonymous + 
                        scholars + 
                        opposition*trust_parliament_pretreatment + 
                        veteran +
                        Dem_Age + Dem_Gender_Others + Dem_Gender_Female +
                        SES_Income + SES_Education + SES_Interest,
                      se_type = "stata", data = d)

model3.4 <- lm_robust(outcome_dispatch ~ anonymous + 
                        scholars + 
                        opposition + 
                        veteran*trust_SDF_pretreatment +
                        Dem_Age + Dem_Gender_Others + Dem_Gender_Female +
                        SES_Income + SES_Education + SES_Interest,
                      se_type = "stata", data = d)

model3.5 <- lm_robust(outcome_dispatch ~ anonymous + 
                        scholars*trust_university_pretreatment + 
                        opposition + 
                        veteran +
                        Dem_Age + Dem_Gender_Others + Dem_Gender_Female +
                        SES_Income + SES_Education + SES_Interest,
                      se_type = "stata", data = d)


cplot(model3.2, x = "trust_scholars_pretreatment", dx = "scholars", what = "effect", rug = FALSE, ylim = c(-0.65, 0.55), 
      main = "Criticism from Constitutional Scholars", xlab = "Confidence in Constitutional Scholars",ylab = "Marginal Effect",
      cex.main =1.5, cex.lab =1.5, cex.lab = 1.5, draw = TRUE)
abline(h=0, lty=2)

cplot(model3.3, x = "trust_parliament_pretreatment", dx = "opposition", what = "effect", rug = FALSE, ylim = c(-0.65, 0.55),
      main = "Criticism from Opposition Party", xlab = "Confidence in the Diet", ylab = "Marginal Effect",
      cex.main =1.5, cex.lab =1.5, cex.lab = 1.5,draw = TRUE)
abline(h=0, lty=2)

cplot(model3.4, x = "trust_SDF_pretreatment", dx = "veteran", what = "effect", rug = FALSE, ylim = c(-0.65, 0.55),
      main = "Criticism from Retired General Officers",xlab = "Confidence in SDF", ylab = "Marginal Effect",
      cex.main =1.5, cex.lab =1.5, cex.lab = 1.5,draw = TRUE)
abline(h=0, lty=2)

cplot(model3.5, x = "trust_university_pretreatment", dx = "scholars", what = "effect", rug = FALSE, ylim = c(-0.65, 0.55),
      main = "Criticism from Constitutional Scholars", xlab = "Confidence in Universities",ylab = "Marginal Effect",
      cex.main =1.5, cex.lab =1.5, cex.lab = 1.5,draw = TRUE)
abline(h=0, lty=2)
```




#Balance check


```{r}
#remotes::install_github("JaehyunSong/BalanceR")  

library(BalanceR)
```


```{r}
d$SES_PoliticalInterest <- d$SES_Interest

d$group <- ifelse(d$control == 1, "Control",
                 ifelse(d$anonymous == 1, "Anonymous",
                 ifelse(d$scholars == 1, "ConstitutionalScholars",
                 ifelse(d$opposition == 1, "Opposition", "FormerGeneral"))))

BlcChk <- BalanceR(data = d, group = group,
                   cov = c(Dem_Age, Dem_Gender_Female, Dem_Gender_Others,
                           SES_Income_300under, SES_Income_300_350, SES_Income_350_400, SES_Income_400_450, SES_Income_450_500,
                           SES_Income_500_550, SES_Income_550_600, SES_Income_600_650, SES_Income_650_700, SES_Income_700_750,
                           SES_Income_750_800, SES_Income_800_850, SES_Income_850_900, SES_Income_900_950, SES_Income_950_1000, SES_Income_1000over,
                           SES_PoliticalInterest,
                           SES_Education_HighSchool, SES_Education_VocationalHighSchool, SES_Education_University, SES_Education_GraduateSchool,
                           ideology,
                           trust_scholars_pretreatment, trust_parliament_pretreatment, trust_SDF_pretreatment))

DF <- print(BlcChk)
```



Visualize the maximum absolute values of standardized biases.
```{r}

plot(BlcChk, simplify = TRUE)

```





